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Abstract 

We prove a completeness result for a class of polynomial solutions of the wave equation called 
wave polynomials and construct generalized wave polynomials, solutions of the Klein-Gordon 
equation with a variable coefficient. Using the transmutation (transformation) operators and 
their recently discovered mapping properties we prove the completeness of the generalized wave 
polynomials and use them for an explicit construction of the solution of the Cauchy problem 
for the Klein-Gordon equation. Based on this result we develop a numerical method for solving 
the Cauchy problem and test its performance. 



1 Introduction 

Let fi C C be a simply connected domain. Due to the Runge approximation theorem any harmonic 
in Q function can be approximated uniformly on any compact subset inside Q by harmonic poly- 
nomials. The harmonic polynomials are linear combinations of the polynomials Re(z — zo) n and 
lm(z — zo) n , re = 0, 1, ... , where zo is an arbitrary point in Q and z is a complex variable. This fact 
reflects the completeness of the system of harmonic polynomials {Re(z — zo) n , lm(z — zo) n }^L in 
the space of all harmonic functions in Q in the sense of the normal convergence. 
Instead of the Laplace equation let us consider the wave equation 

w xx -w u = (1) 
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supported by DFFD, Ukraine (GP/F32/030) and by SNSF, Switzerland (JRP IZ73Z0 of SCOPES 2009-2012). The 
research of Sebastien Tremblay is partly supported by grant from NSERC of Canada. 



and instead of the complex imaginary unit let us introduce the hyperbolic imaginary unit: j 2 = 1. 
Let z denote the hyperbolic variable z = x + jt [22], [28]. Analogously to the elliptic case the 
system of polynomials 

{Re(x + jt) n and Im(x + jt) n }™ =0 (2) 

is an infinite system of solutions of the wave equation. Up to now, to our best knowledge, no 
corresponding completeness result has been obtained. We call the polynomials ^ and their finite 
linear combinations wave polynomials, and one of the first results of the present work is a Runge- 
type theorem establishing that any regular solution of |I| in a closed square R with the vertices 
(±26, 0) and (0, ±26), 6 > can be uniformly approximated on R by the wave polynomials. This 
theorem is auxiliary for obtaining a similar result for solutions of the Klein-Gordon equation with 
a variable coefficient 

Uxx ~ utt ~ q(x)u = (3) 

which we consider next. The construction of an infinite system of solutions similar to the wave 
polynomials was done in |19j with the aid of L. Bers' results on pseudoanalytic formal powers 
[2] extended onto the hyperbolic situation. Similarly to the wave polynomials these generalized 
wave polynomials are components of formal powers, solutions of a corresponding hyperbolic 
Vekua equation which locally behave as powers of z = x + jt but in general are not of course 
powers. Using recent results from [3j on mapping properties of transmutation operators we show 
that the generalized wave polynomials are images of the wave polynomials under the action of a 
transmutation operator. Due to the uniform boundedness of the transmutation operator and of 
its inverse several useful properties of the wave polynomials are preserved also in the case of their 
generalizations. In particular, the expansion theorem and the Runge-type theorem result to be 
valid. 

All these observations lead to a new representation for the solution of the Cauchy problem for 
Q. It is based on the expansion of the Cauchy data into series in terms of a certain system of 
functions {^fcl^Lg wn i cn are introduced as recursive integrals and arise in the spectral parameter 
power series (SPPS) representation for solutions of Sturm-Liouville equations [14] . [18] . In [16] 
a completeness of {<fk}fc=o in L2 was proved. In [T7] this result was obtained for the space of 
continuous and piecewise continuously differentiable functions. Here we show that the completeness 
of {</?fc}fcLo in the space of continuous functions directly follows from the mapping properties of 
the transmutation operator and the Weierstrass approximation theorem. In [T7] it was shown that 
several classical results from the theory of power series can be generalized onto the series in terms 
of the functions (fk, including the Taylor formula. Here we present several new results on the 
approximation of continuous functions by linear combinations of functions (fk- In particular, we 
show that the system of functions {</?fc}fcLo i n a real-valued case is a Tchebyshev system, prove a 
direct and an inverse approximation theorems and study algorithms for such approximation. 

Using the results on the approximation by functions tpk we propose a numerical method for 
solving the Cauchy problem for ^ and illustrate its performance with several test examples. Once 
the Cauchy data are approximated by functions (pj~, the approximate solution of the Cauchy problem 
is written in a closed form. As for t > the approximate solution is an exact solution of equation 
Q the only task consists in a good approximation of the Cauchy data. We show that in fact with 
relatively few functions <p>k involved, a remarkable accuracy is achieved. 
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2 Wave polynomials 

Let us consider the wave equation 

d 2 d 2 

□» = 0, V: = — 2 - W2 (4) 
and the following infinite family of its polynomial solutions 

{Re(x + jt) n and Im(x + jt) n }™ =0 (5) 

where j is a hyperbolic imaginary unit, j 2 = 1. 
It is easy to see that 

Re(x + jt) n = - ((x + t) n + (x - t) n ) and Im(x + jt) n = ~ ((z + t) n - (x - t)") . (6) 
Let us reorder these polynomials as follows 

po(x,t) = l, pi(x, t) = Re(x + jt) = x, p2(x,t) = lm(x + jt) = t, 
p 3 (x, t) = Re(x + jt) 2 = x 2 + t 2 , p 4 (x, t) = Im(x + jt) 2 = 2xt,. . . . 

The obtained family of solutions of Q will be called wave polynomials. It is convenient to write 
them also in the following form 

( m + l 

2 /m+l 



p (x,t) = l, p m (x,t) 



2 / m+l \ m+1 

^ v fe / ^ ^ m odd ' 



even fc=0 

m 

2 /m 



(7) 



( 2 jx™ m even. 



V odd fc=l 

Consider equation Q together with the following Goursat conditions 

u> = <p(x) for x — i = and u; = ip(x) for x + t = (—6 < x < 6), 

assuming additionally that (f(0) = ip(0). It is well known (see, e.g., [3Q|. 4.1.1-9.]) that for <p and ifi 
belonging to C 2 [— 6, b] the solution of the Goursat problem exists, is unique and has the form 

w (x,t) = V (i±*)+^(^)- V »(0). (8) 

Its domain of definition is a closed square R with the vertices (±26, 0) and (0, ±26). 

Proposition 1 Let the boundary data (p and ip be real- analytic functions with the corresponding 
power series expansions 

oo oo 

ip(x) = ^a„x n and ijj(x) = ^/3 n x n , (9) 

n=0 n=0 

uniformly convergent on [—6,6] and satisfying necessary condition yj(0) = ip(0), i.e., ceo = Po- Then 
the unique solution of the Goursat problem has the form 

X,t) = a p {x,t) + 2^ ^ P2n-l{X,t) + 2^ ^ P2n(X,t) 

n=l n=l 

where the series converge uniformly in R. 
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Proof. From ^ we have 



P2n-i(x,t) = ) 2 {{x + t) n + (x-t) n ) and pan(x,t) = \{{x + t) n - (x - t) n ) 

and hence 

(x + t)™ = p 2n -l(x,t) +P2n{x,t) and (x-t)™ =p 2 n-l(^,t) ~ P2n{x,t), 71=1,2,... 

From (Tsl) we obtain that the solution of the Goursat problem has the form 



• (10) 



w(x,t) = a + ay 



(re + i) r 



(s - t) n 



n=l n=l 

Substitution of the relations (flO|) gives us the equalities 



w[x 



. p 2n -l(x,t) +p 2n (x,t) ^ 

,*) = a + }^a n + ^ 



n=l 



n=l 



ggn-lCgjj) -P2n(x,t) 

2 n 



a Po(s,i) + 2^ ™ P2n-l(X,t) + 2^ ^ P2n{X,t). 



n=l 



n=l 



Remark 2 From this proposition we obtain that the wave polynomials represent a basis in the 
linear space of solutions of the wave equation which admit a uniformly convergent in R power 
series expansion with the center in the origin. Indeed, consider any such solution of Q in R. Rs 
values on the lines x — t = and x + 1 = admit uniformly convergent power series expansion of 
the form ^ . According to the proposition the considered solution can be represented as a uniformly 
convergent series of the wave polynomials. 

Let us prove the completeness of the wave polynomials in the linear space of regular solutions 
of the wave equation with respect to the maximum norm. 

Theorem 3 Let w G C 2 (R) be a solution of the wave equation Q in R. Then there exists a 
sequence of wave polynomials Pn = Yln=o a nPn uniformly convergent to w in R. 

Proof. We need to prove that for any e > there exist such a number ./V and coefficients 
a n , n = 0,1, ...N that \w(x, t) — Pjy(x, t)\ < e for any point (x,t) £ R. Let w = f(x) for 
x — t = and w = ip(x) for x + t = (—b < x < b). We choose e > and such > 
that e = 2e\ + 2^2- According to the Weierstrass theorem there exists such number N and such 
polynomials p\ and P2 of order not greater than N that \<p(x) — pi(x)\ < E\ and \ip(x) — P2(x)\ < £2 
(—6 < x < b). We consider polynomials q\{x) = pi(x)—p\(0)+ip[0) and q^ix) = P2(x)—p2(0)+ip{0) 
satisfying the condition q\(0) = 92(0) = <p(0). Due to Proposition [l] the unique solution w of the 
Goursat problem with the boundary data q\ and q2 has the form w = Pj\[(x,t) where Pw(x,t) = 
qi (^) +Q2 (^) -9i(0). Consider 



\w(x, t) — w(x, t)\ 



< 



< 



+ 



w(x, t) - 
x + t 
2 

x + t 



PN{x,t)\ 

X + t 



1> 



2 

x — t 



pi 



P2 



2 

X + t 

2 

x — t 



+ 



X 



qi 



x — t 
2 



+ 1^(0) -pi (0)| 

+ 1-0(0) -P2(0)| < 2ei + 2e 2 = e. 
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3 Transmutation operators and their action on powers of the in- 
dependent variable 

3.1 Systems of recursive integrals 

Let / G C 2 (a,6) n C^a, 6] be a complex valued function and f(x) / for any x G [a, b]. The 
interval (a, b) is supposed to be finite. Let us consider the following auxiliary functions 

X<®(x) = X(°\x) = 1, (11) 
xW(x) = n f X X^(s) {f 2 (s)) { - ir " ds, (12) 

X^(x) = n f X^Xs) {f 2 {s)) { - ir ds, (13) 

where xq is an arbitrary fixed point in [a, b]. We introduce the infinite system of functions {^fcl^Lg 
defined as follows 

ff(x)X( k \x), kodd, 
Vk{x) = < ~ (k) (14) 

A: even, 



where the definition of X^ and X^ is given by ( 11 )— ( 13 ) with xo being an arbitrary point of the 
interval [a, b]. 

Example 4 Let / = 1, a = 0, 6 = 1. T/ien it is easy to see that choosing xq = we have 
'fk(x) = x k , k G No where by No we denote the set of non-negative integers. 

In [16] it was shown that the system {^fcjfcLo is complete in L,2(a, b) and in [T7] its completeness 
in the space of continuous and piecewise continuously differentiable functions with respect to the 
maximum norm was obtained and the corresponding series expansions in terms of the functions (pk 



were studied. The completeness in the space C[a, 6] is shown in the Proposition 27 



The system ( 14 ) is closely related to the notion of the L-basis introduced and studied in [8] . Here 
the letter L corresponds to a linear ordinary differential operator. This becomes more transparent 
from the following result obtained in [14J (for additional details and simpler proof see [15] and [18J) 
establishing the relation of the system of functions {^fcjfcLo to Sturm-Liouville equations. 

Theorem 5 (|14|) Let q be a continuous complex valued function of an independent real variable 
x £ [a, 6], A be an arbitrary complex number. Suppose there exists a solution f of the equation 

f" ~qf = (15) 

on (a, 6) such that f G C 2 [a, b] and f ^ on [a, 6]. Then the general solution of the equation 

u" — qu = Xu (16) 

on (a, 6) has the form 

U = C\U\ + c 2 u 2 
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where c\ and c 2 are arbitrary complex constants, 



A 

Ui 



Ep)T^ and ^ = ^2^—^^ (17) 



and both series converge uniformly on [a, b] . 

The solutions u\ and U2 satisfy the initial conditions 

ui{x ) = f(x ), u[(x ) = f'(x ), 

u 2 (x ) = 0, u' 2 (x ) = l/f(x ). 

Together with the family of functions {i^fcjfcLo we consider a dual system of recursive integrals 



defined by the following relations involving the "second half of the formal powers ( 11 )— ( 13 ) 

k odd, 



tpk(x) 



' X( k \x) 



(18) 

k even. 



3.2 Generalized derivatives and generalized Taylor series 

In |17| a notion of the generalized derivative was introduced which alows one to extend the theory 
of power series onto the series in terms of the functions ipi~ (the formal power series). Here we 
slightly modify the definition introduced in [17J. This modification simplifies formulas involving 
the generalized derivatives and reflects a better understanding of the nature of the functions <pk 
and ipk in the light of application of transmutation operators. We assume that the complex- valued 
function / is continuous on [a, b], f(x) ^ for any x G [a, b] and f{xo) = 1. 

Definition 6 The generalized derivatives or the f -derivatives of a function g are defined by the 
following relations whenever they make sense. The generalized derivative of a zero order coincides 
with the function g, d^[g](x) = g(x). The generalized derivatives of higher orders are defined as 

follows 4[g] = f^ h - 1 i; (f(-V k 4_ 1 [ g }), 
That is, 

4[g] 



A; = 1,2, 



f£{H-ii9]h k odd, 



even. 



Remark 7 Let f be a solution of (15) satisfying the conditions of Theorem^ Then the correspond- 
ing differential operator can be factorized in the following way L = — q(x) = j-^ (j 2 -^.j • J ■ 
This factorization sometimes is called the Poly a factorization (see \1 1^). We see from it that L = d 2 - 
The generalized derivative d{ = (^j coincides with the Darboux transformation (see, e.g., 

Remark 8 It is easy to see that 

d{ifk = kij) k -i, A; = 1,2,..., 
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d{ipk = k(k— l)<Pk-2, fe = 2,3, — 

and 

f f f 

d{(po = d 2 (fo = d 2 (fi = 0. 

Remark 9 Consideration of the 1/ /-derivatives defined according to Definition^leads to the dual 
relations 

l it 

d\ i>k = k<Pk-i, A; = 1,2,... 

and 

i /f d ( 1 d \ d 2 

d * =f Tx\ y pd-J-) = d^- qD{x) 

f s2 



where the potential qu is a superpartner of q defined by the equality qr> = —q + 2[h-\ (see 



Subsection 3.3) 



Definition 10 Functions of the form 

n 

p n{x) = ^a k tp k {x) (19) 

k=0 

where a k , k = 0, 1, . . . , n are complex numbers are called f -polynomials of the order n. 

In a complete similarity to the fact that the coefficients of a polynomial X^fc=o a k( x ~ x o) k can 
be expressed through its value and the values of its derivatives at the point xq, the coefficients of an 
/-polynomial are determined by the values of P n and of its /-derivatives at xq (at the initial point 



of integration in (12), (13)). Indeed, a simple calculation using Remark [8] gives us the following 
result 

4[Pn](x ) 

ak = — a — • 

Let us consider a function g possessing at the point xq the /-derivatives of all orders up to the 
order n. More precisely this means that the function g is defined and possesses the /-derivatives 
of all orders up to the order re — 1 in some segment [a, b] containing the point xq and additionally 
there exists the n-th /-derivative of g at the point xq. In relation with the function g, we introduce 



an /-polynomial of the form (19) with the coefficients 



d{[0](a?o) 
ak = ^— 

According to the previous observation, this /-polynomial together with its /-derivatives at xq up 
to the order re take the same values as the function g and its respective /-derivatives, d^ k [P n ]{xo) = 
d{[g](xo), k = 0, 1, . . . , n. We are interested in estimating the difference between P n (x) and g(x) 
for x ^ Xq. 

Theorem 11 (Generalized Taylor theorem with the Peano form of the remainder term) 

Let the function g possesses at the point xq the f -derivatives of all orders up to the order n and f 
be a continuously differentiable function in a neighborhood of xq. Then 

9( X ) = Z^ u M x ) + o{{x-x ) ). 

k=0 
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Proof. Consider the difference r(x) = g{x) — P n (x). We have 

r(xo) = d{[r](x ) = ■■■ = d f n [r](x ) = 0. (20) 



Let us prove by induction that if a function r satisfies the conditions (20) or the conditions 

r(x ) = d{ /f [r}(x ) = ■■■ = d^ f [r](x ) = 0, (21) 

then necessarily r(x) = o((x — xo) n ). 

For n = 1 this assertion has the form: if the function r(x) possessing at xq the first /-derivative 

fulfills the conditions r(xo) = d[[r](xo) = or possessing at xq the first l//-derivative fulfills the 

1/ f 

conditions r(xo) = <i 1 [r](sco) = then r(x) = o[x — xo). Its validity can be verified directly. In 
the first case we have 

( \ I Mft \ r(x)_ _ r(x ) 

iim JM. = f{xo) lim r M!M = /(xo) lim m jm = o 

x-^xq x — Xo x-+xq x — Xo x->x X — Xo 

due to the condition d{[r](xo) = 0, and in the second case the proof is completely similar. 



Assume that the assertion is true for some n > 1. Due to the symmetry of (20) and (21) 
it is enough to prove that if for a function r(x) possessing at xo the /-derivatives up to the 
order n + 1 the following conditions are fulfilled r(xo) = d{[r](xo) = ■ ■ ■ = d{ l+1 [r](xo) = then 
r(x) = o((x — xo) n+1 )- For this we observe that r(x) fulfills the conditions (20) meanwhile d{[r] 
fulfills the conditions ( |2l[ ) and hence by the assumption we have r(x) = o((x — xo) n ) and d[[r](x) = 
o{ix — xo) n ). Notice that by the mean value theorem 

r(x) =r(x) — r(xo) = (Rer'(ci) + zlmr'(c2)) (x — xq) 

Re(d{[r]( Cl ) + ^r( Cl )) +il m (d{[r](c 2 ) + ^r(c 2 ))) (x - x ), 

where ci and c 2 are located between xq and x. As \ci t 2 — xq\ < \x — xq\, then d{[r\{ci t 2) = o((ci ) 2 — 
xo) n ) = o((x — xo) n ) and r{c\p) = o((ci j2 — xo) n ) = o((x — »o) n )- Thus, we obtain r(x) = 
o((x-x ) n+1 ). ■ 

Under an additional condition that / is real valued we obtain the following result by applying 
the reasoning from |17| . 

Theorem 12 (Generalized Taylor theorem with the Lagrange form of the remainder) 

Let the real-valued function g possesses on the segment [xo,b] continuous f -derivatives of all orders 
up to the order n and there exists a bounded (n + l)-th f -derivative of g on (xo,b). Let f be a 
real-valued, continuously differentiate function in [xo,6]. Then for any x G [xo,6] there exists a 
number c between xq and x such that 

( v ^ d{(g){xo) d f n+1 (g){c) 

9^ = 2^ jfcj Vkix) + (n+1) , <Pn+l(x). 

Proof. The proof is a simple adaptation of the proof from [17] according to the modified definition 
of generalized derivatives. All the steps and reasonings do not essentially change. ■ 

Obviously, the classical Taylor theorem with the Lagrange form of the remainder term is a 
special case of theorem |12| when / = 1. 
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3.3 Transmutation operators 

We give a definition of a transmutation operator from [21] which is a modification of the definition 
given by Levitan [23] adapted to the purposes of the present work. Let E be a linear topological 
space and E\ its linear subspace (not necessarily closed). Let A and B be linear operators: E% — > E. 

Definition 13 A linear invertible operator T defined on the whole E such that E\ is invariant 
under the action of T is called a transmutation operator for the pair of operators A and B if it 
fulfills the following two conditions. 

1. Both the operator T and its inverse T~ l are continuous in E; 

2. The following operator equality is valid 

AT = TB (22) 

or which is the same 

A = TBT~ l . 



Very often in literature the transmutation operators are called the transformation operators. 
Here we keep the original term introduced by Delsarte and Lions [5J. Our main interest concerns 
the situation when A = — 4-% + q(x), B = — -jiz, and q is a continuous complex- valued function. 
Hence for our purposes it will be sufficient to consider the functional space E = C[a, b] with the 
topology of uniform convergence and its subspace E\ consisting of functions from C 2 [a,b). One 
of the possibilities to introduce a transmutation operator on E was considered by Lions [21] and 
later on in other references (see, e.g., [25]), and consists in constructing a Volterra integral operator 
corresponding to a midpoint of the segment of interest. As we begin with this transmutation 
operator it is convenient to consider a symmetric segment [—b, b] and hence the functional space 
E = C[—b,b]. It is worth mentioning that other well known ways to construct the transmutation 
operators (see, e.g., [23], [38]) imply imposing initial conditions on the functions and consequently 
lead to transmutation operators satisfying (22) only on subclasses of E\. 



Thus, we consider the space E = C[—b,b] and an operator of transmutation for the defined 
above A and B can be realized in the form (see, e.g., [23] and [25]) of a Volterra integral operator 



Tu{x) = u{x)+ K(x,t)u(t)dt (23) 

J — X 

where K(x,t) is a unique solution of the Goursat problem 

f d 2 \ d 2 

{d^-^ x )) K ( x ^ = W Kix ^ (24) 

i r 

K(x,x) = - q(s)ds, K(x,-x)=0. (25) 



o 



2 

In [3] the following mapping properties of the operator T were proved 



Theorem 14 (|3j) Let q be a continuous complex valued function of an independent real variable 
x £ [—b,b] for which there exists a particular solution f of (15) such that f E C 2 [— b, b], f ^ 



on [—b,b] and normalized as /(0) = 1. Denote h := /'(0) G C. Suppose T is the operator defined 
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by (23) where the kernel K is a solution of the problem (24), (25) and (p^, k G No are functions 
defined by (14). Then the following equalities hold 



and 



h 

k + 1 



T\x k 



ifk+i = T[x k 



k odd 



k even. 



Remark 15 Let f be the solution of (15) satisfying the initial conditions 

/(0) = 1 and /'(0) = 0. 



(26) 
(27) 

(28) 



If it does not vanish on [—b,b] then from Theorem 14 we obtain that tpk = T[x ] for any k G No- 
In general, of course there is no guaranty that the solution with such initial values would have no 
zero on [—6,6], and hence the operator T transmutes the powers of x into <pk whose construction is 
based on the solution f satisfying (28) only in some neighborhood of the origin. 



In [3] it was shown that given a system of functions {(^fcj^Q defined by (14) where / is any 
particular solution of (l5) such that / G C 2 [-6,6], / / on [-6,6] and /(0) = 1, /'(0) = h G C, 
the operator T can be modified in such a way that the new transmutation operator will map x k to 
(fk for any k G Nq. 



Theorem 16 ([3], [20]) Under the conditions of Theorem 14 the operator 



Tu(x) = u(x) + 

with the kernel defined by 

K(x,t;h) 

transforms x k into (fk{x) for any k G No and 

d 2 



h s h 

- + K(x,t) + - 



K(x,t;h)u(t)dt 



(K(x,s) - K(x,-s))ds 



dx 2 



+ q(x) T[u] 



dx 2 



(«) 



(29) 



(30) 



(31) 



for any u G C 2 [— 6, 6]. 

Moreover, if the potential q G C l [— 6, 6], then the kernel K.(x,t;h) is a unique solution of the 
Goursat problem 

92 ^ 82 ,32) 

:. (33) 



K dx 2 
K(x, x; h) = 



q{x) ) K(x,t;h) = -^K(x,t;h), 



h 1 

2 + 2 



q(s)ds, 



K(x, —x; h) 



This theorem was proved in [3] under an additional assumption that the potential q must be 
continuously differentiable, and in [20] it was shown that this assumption was superfluous due to 
new observed relations (38), (39) given below between the transmutations for Darboux transformed 
Schrddinger operators. The last statement of this theorem was proved in [20], also the interested 
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reader may find in |20j necessary changes regarding the case when q S C[—b,b]. For brevity, we 
omit these details in the present article. 

In the following sections we use both the transmutation operator T and its inverse T _1 , and 
the norms of these operators appear in many estimates. Hence it is natural to obtain convenient 
estimates for the norms. Remind that in [20] it was mentioned that to define the transmutation 
operator T, we need to know its integral kernel in the domain < |i| < |x| < 6. But the Goursat 



problem (32)-(33) is also well-posed and allows to determine the kernel K.(x,t;h) in the domain 
< \x\ < \t\ < b. Thus we may assume that the integral kernel K(x, t; h) is known in the square 
|x| < 6, |t| < 6. In such case, there is a simple representation of the inverse operator T . 

Theorem 17 (|20j) The inverse operator T _1 can be represented as the Volterra integral operator 

T~ 1 u(x) = u(x) — [ K(t,x;h)u(t)dt. (34) 



Both T and T _1 are obviously bounded as operators from C[—b,b] to itself. The estimates 
for their norms depend on the estimates for the integral kernels, e.g., for ||T|| we have ||T|| < 
1 + 26 max |K(x, t; h)\. Some estimates for the integral kernel K(x,t) can be found in [25] . From 



them corresponding estimates for the kernel K(x,i;/i) can be obtained using (30). However, the 
growth with the increase of the interval of mentioned estimates is immensely fast even for the 
simplest potentials. We adapt the general method of successive approximations for solving Goursat 
problems (see, e.g. [39]) to obtain better estimates for the kernel K.(x,t;h). 

Proposition 18 Let q be a continuous complex valued function of an independent real variable 
x G [—6,6]. Then the kernel JZ(x,t;h) in the square \x\ < b, \t\ < 6 satisfies the following estimate 



\h\ 



\K(x,t;h)\<!fl (J. 



cr 



' c 2r 



t 2 \hy C \x 2 -t 2 \) 



(35) 



where c = max[_j, ^ \q{x)\ and Iq and I\ are modified Bessel functions of the first kind. 



Remark 19 Note that for the case of operator d 2 — c with constant potential c > and h > the 
exact kernel of transmutation operator in the domain < t < x < 6 coincides with the right-hand 



side of (35), see |2|/. 



Proof. The proof follows the proof from [3S]. First, we introduce the function H(u,v) := K(u + 
v, u — v; h). It satisfies the Goursat problem (see [20J) 

d 2 H{u,v) 



H(u,0) 



h 1 

2 + 2 



q(u + v)H(u,v), 
q(s)ds, H(0,v) 



du dv 

h 

Jo 2 
in the domain \u\ + \v\ < b. It is worth mentioning that despite the kernel K(x,£;/i) is not the 



(36) 
(37) 



classical solution of the problem (32)-(33) in the case when q S C[— 6, 6], nevertheless the function 
H(u, v) is a classical solution of the problem (36)-(37), see [20]. Define G:=^. Then the Goursat 



problem (|36h— ( 37 ) is equivalent to the system of integral equations 



'H{u,v) = | + J " G{u\v)dv! 
G(u, v) = \q{u) + /* q(u + v')H(u, v') dv' . 
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Applying the successive approximations method for this system, we obtain 



\H(u,v)\ < 



\h\ 



EC \UV\ 1 



c k+1 \u\ k+1 \v\ k 



2 ^ k\k\ 

k=0 



2f ( k + ^ 



which coincides with (35). 



Since the function I\(x)/x is monotone increasing for x > 0, we obtain 



^c\x 2 -t 2 \h ( ^c\x 2 -t 2 \) 



c\x + t 



''<^< 26c W^) 

t 2 \ - byfc 



CI' 



for |x| < 6 and \t\ < 6, and the following estimate for the norms of transmutation operator and of 



its inverse immediately follows from Proposition 18 
Corollary 20 The following estimate holds 



maxdlTHJT- 1 !!} < 1 + b (\h\I (b^c) + 2y/Blx(by/c)) , 
where c = ma.x^_ b b ^ \q(x)\ and Iq and I\ are modified Bessel functions of the first kind. 

Together with the operator —q(x) let us consider a Darboux associated operator 



dx- 



■Qd(x) 



with the potential defined by the equality qr> 



-g + 2('=7 



where / G C 2 [— 6, b] is a solution 



of (15), / ^ on [—6,6], /(0) = 1 and h = /'(0) G C. In [20J explicit formulas were obtained 
for the kernel Kjj(x,t; —h) in terms of K(x, t; h), where K£>(x,i; —h) is the integral kernel of the 
transmutation operator Try which satisfies the equality 



dx 2 



+ q D {x) T D [u] =T D 



dx 2 



u 



for any u G C 2 [— 6, 6] and transforms x k into the functions ipk{x), k 6 No defined by the relations 
(18). Note that ipo is obviously a solution of ^— + Vo = with the initial values "00(0) = 1 

.%(Q) = -h. 

The operator has the form |20j 



and 



with the kernel 



T D [u](x) = u(x) + / K D (x,t;-h)u(t)dt, 

J —X 

K D (x, t; -h) = --^ ( J* 8 t K(s, t; h)f(s) ds + ^f(-t) 



and the following operator equalities hold on C 1 [— 6, 6]: 

dx dx 



dx f 



-T — 

/ ° dx' 



(38) 
(39) 



These commutation equalities involving the operators of transmutation and derivatives together 
with the property of the transmutation operators that if u G C 1 [— 6, 6] then T _1 u G C 1 [— 6, 6], see 
Theorem 17 lead to the following useful statement. 
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Proposition 21 ( |21j ) Let u G C n [— b, b] and g = Tu. Then there exist the first n f -derivatives 
of g on [—b, b], and the following equalities hold for < k < n 

d{{g) = T D u {k \ k odd, 

and 

d{(g) = Tu^, k even. 

The inverse statement, i.e., if there exist the first n f -derivatives of g on [—b, b], then u = T -1 ^ G 
C n [— b, b] is also true. 



4 Generalized wave polynomials 

Let us consider the following Klein-Gordon equation with a position dependent mass 

d 2 



(dx 2 



(x) u(x,t) 



dt 2 



u(x, t) 



(40) 



where we assume that q : [—b,b] — > C and q G C[— b, b]. Suppose there exists a particular solution 
/ of equation (15) such that / G C 2 [— b, b] and / ^ on [—b,b\. We normalize it as /(0) = 1 and 
set h := /'(0). 

Consider the system of functions {<^fc}^Q_defined by (14) with xo = 0. Then due to Theorem 
(fk{x) = Tx fc for any k G No and due to Q7J) we obtain that the functions 
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u = f(x), u m (x,t) 



m+1 

2 

E 

even k=0 

m 
2 

E 

^ odd k=l 



m+1 s 
2 („\+k 

k 



(f m+i k (x)t k , m odd, 



(41) 



2 )^^_ k {x)t k , 



rn even, 



are solutions of (40) for any —b < x < b and — oo < t < oo. Indeed, we have that 

Tp m for every m G Nq. 



U r , 



(42) 



Moreover, the functions u m arise also as scalar (real, when / is real valued) parts of hyperbolic 
pseudoanalytic formal powers corresponding to the generating pair (f,j/f) where j is a hyperbolic 
imaginary unit, j 2 = 1 (see [H], [T5"]). 

Equalities ( 42 ) together with the completeness of the wave polynomials (Theorem [3]) and the 
boundedness of T and T" 1 imply the completeness of the generalized wave polynomials u m in 
the linear space of regular solutions of (40). 

Theorem 22 Let u G C 2 (R) be a solution of (40) in R where R is a square with the vertices 
(±6, 0) and (0, ±6). Then there exists a sequence of generalized wave polynomials Un = X^=o a nU n 
uniformly convergent to u in R. 

Proof. We have that u = Tw where w is a C 2 -solution of Q and due to Theorem [3] for any ei > 
there exists a wave polynomial P/y such that rnax^ \w — P/v| < £\. Thus, max-^ |u — TP/v| = 
max-jj \Tw — TPj\r| < E\C = e. Here the constant C depends only on the kernel K(x, t; h). ■ 
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Remark 23 When t = the following relations are valid 



and 



u m (x,0) 

du m (x,0) 
Ft " 



(f m+i (x), m odd 



0, 



m even 



0, m odd 

(pni-i(x), m even. 



in 

2 ' 2 



These relations follow directly from the definition (41). One can write them also as follows 



u 2n -i(x, 0) = ip n (x), u 2n (x, 0) = 0, 



anc 



du 2n -i{x,0) 



0. 



du 2n (x,0) 



For uq we have 



dt ' dt 

Uq(x, 0) = f(x) = ipo(x) and 



n<p n -i(x), for n = 1,2, ... . 
duo(x, 0) 



dt 



0. 



5 Solution of the Cauchy problem 

Consider the following initial value problem 

Du-q(x)u = Q, -b<x<b, t>0 

u(x,0) = g(x), u t (x,0) = h(x) 

which for q £ C[—b,b], g G C 2 [—b,b], h G C 1 [— 6, 6] possesses a unique 
solution (see, e.g., [391 Sect. 15.4]) in the triangle with the vertices 
(±6,0) and (0,6) (see illustration). For the convenience, later in this 
article we denote this triangle by the symbol ▲. We assume that q 
satisfies the conditions of Theorem 14 and begin with the additional 
assumption that the functions g and h admit uniformly convergent 
series expansions in terms of the functions tfik, 



9{x) = ^2 "Wkix) and h(x) = ^ Pk¥k{x) 



k=0 



k=0 



We look for a solution of the problem ( 43 ) , ( 44 ) in the form 



U[X 



t) = y2 a n u n(x,t). 



n=0 



Then we have (see Remark 23 ) 



u(x, 0) = a ip (x) + ^2 a 2n -lU2n-l(.x, 0) = a ip (x) + ^ a 2n _i</?„(3 



(43) 
(44) 





t' 


b 












-6 






b x 


x). 






(45) 



(46) 



n=l 



n=l 
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and 



u t 



(x,0) = y^a2n 



n=l 



9n 2 n(^,0) 

at 



^a 2 „n( ( c) n _i(j;) = ^a 2 (fc+i)(^+ l)Vfc(a). 



n=l 



fc=0 



Thus, if a solution of the problem (43), (44) in the form (46) exists, the expansion coefficients are 



obtained directly from the coefficients in the expansions (45) as follows 



ao = «0> «2n-l 



Otr, 



n 



1,2,... 



and 



a 2(n+l) 



fin 
71+1' 



n 



0,1,2, 



(47) 



The following natural questions arise. Under which conditions given functions g and h are 



representable in the form (45) and whether such series expansion is unique? Can one guarantee 



the uniform convergence of the series ( 46 ) and that of its first and second derivatives in a domain 



of interest? In what follows we address these questions and show that the described scheme also 



leads to a powerful numerical technique for solving the initial value problems for equation ( 43 ) 



Proposition 24 A continuous complex-valued function g defined on [—6, 6] admits a series ex- 
pansion of the form g(x) = J2h=o a k ( Pk(x) uniformly convergent on [—b,b] if and only if there 
exists a complex-valued function g defined on [— b, b] such that g = Tg, g(x) = Y^k=o a k xk an< ^ ^ e 
power series converges uniformly on [—6,6]. The expansion coefficients are uniquely defined by the 
equalities 



4(g)(0) gW(p) 



k\ 



(48) 



Proof. The proof of the represent ability of g in the form of a uniformly convergent series 
X^fcLo a k i fk(x) follows from the uniform boundedness of the Volterra integral operators T and 
T _1 . The linearity of these integral operators together with the fact that T = ip^ (Theorem 



16) gives us the equality between the coefficients of the corresponding series expansions of g and 
g. The equality d{ (g)(0) = g( k \0), k = 0, 1, . . . is a consequence of Proposition 
observation that at the origin T [u] (0) = Td [ u ] (0) = u(0) for any continuous function u. 



and of the 



Proposition 25 Suppose that g £ C[—b, 6] and for any £ N and x 6 (—6, 6) there exists the 
generalized derivative d k (g)(x) such that for any [—a,a\ C (—6,6) the inequality holds 

4(g)\<C(a;k)^ 

where the constants C(a;k) do not depend on x and the sequence C(a;k) is of a subexponential 
growth dim y/C(a; k) <1). Then on (—6,6) the function g admits a normally convergent general- 
ized Taylor series expansion 



d( x ) = ^Z a Wk(3 



(49) 



fc=0 



and cik = dl(g)(0)/k\. 



Proof. Under the conditions of the proposition consider g = T 1 g. From Proposition 21 we have 
that g G C°°(-6,6) and 

k\ 



(k) 



<MC(a;k) ¥ 
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where M 



= max { | 
< 1 1 T — 1 1 1 



jT^ 1 ]!}. Indeed, considering, e.g., an even k we obtain \g 



(*)| 



di(g) 1 1 'JL" ' || max 
Prom t 
the form g(x j — ^ k 
Taylor series expansion (49) 



4(9) 



< T 



-i| 



C(a; k)H and analogously for an odd k. 



lis we obtain that g admits on (—6, 6) a normally convergent Taylor series expansion of 
k and due to Proposition 



Lfc=0 a k x 



24 



g admits a normally convergent generalized 



Proposition 26 Let the initial data g and h admit uniformly convergent series expansions of the 



form (45) on [—b,b]. Then the unique (classical) solution of the Cauchy problem (43), (44) in ▲ 



(47). 



has the form (46) which is uniformly convergent in ▲. The expansion coefficients are defined by 



Proof. As was previously shown if the series ( 46 ) together with the series corresponding to the 



first and second partial derivatives are normally convergent it satisfies equation ( 43 ) as well as the 



Using ( 42 ) we have 



conditions (44). Thus, it remains to prove the uniform convergence of the involved series. 

< ||T|| • max|p n |, and from ^ we obtain 



u,, 



< b n II T I 



in the triangle ▲. Now, taking into account the uniform convergence of the series (45) on [—6, b] we 
obtain the uniform convergence of the series (46) in ▲. The series corresponding to the first and 



second partial derivatives can be majorized in a similar way with the aid of Remarks [7] and [8j ■ 
As was mentioned above in [T7] it was proved that any continuous and piecewise continuously 
differentiable function on [—b, b] can be approximated arbitrarily closely by a finite linear combi- 
nation of the functions (ff;. The existence of a transmutation operator allows to show that the 
condition of piecewise continuous differentiability is superfluous and provides a simple proof of the 
following proposition. 



Proposition 27 Under the conditions of Theorem 14 the system {^fcj^g is complete in C[—b,b], 
any continuous function on [—b, b] can be approximated arbitrarily closely by a finite linear 



i.e., 



combination of the functions ip^. 

Proof. The proof immediately follows from the existence of the transmutation operator, Theorem 
16 and the Weierstrass approximation theorem. ■ 

Thus, even when it is not possible to guarantee the representability of the initial data g and 



h in the form of uniformly convergent series (45), they can be approximated by corresponding /- 



polynomials. The following statement gives us an estimate of the accuracy of the solution u of the 



problem (43), (44) approximated by a solution un corresponding to the approximated initial data. 



Proposition 28 Let P n (x) = X^fc=o a fc^( x ) an< ^ Qn-i(x) = Sfc=o Pk<Pk(x) be f -polynomials, 
approximating the functions g and h respectively on [—6,6] in such a way that max|<7 — P n \ < e± 
and max | h — Q n -i \ < £2- Let uw(x,t) = Ylk=o a k u k(x, t), N = 2n where ao = olq, a2m-i = %, 
for m = 1,2, ... ,n and «2(m+i) = ^pi; f or m = 0,1, . . . ,n — 1. Then 



max \u 
(x,t)eA 



u N \ < || T || IIT- 1 !! (ei +e 2 6) 



(50) 
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Proof. Notice that ujy is a solution of the Cauchy problem for equation (43) with the initial 
conditions un(x,0) = P n {x), UN,t(x,0) = Q n -i(x), x £ [—b,b]. Consider the functions u = 
T _1 [ii] and un = T _1 [ujv];_ They solve the wave^equation Q and satisfy the initial conditions 
u(x,0) = g(x), ut(x,0) = h(x) and un(x,0) = P n (x), UN,t{x,0) = Q n -i(x), x 6 [— b, b] where 
the tilde indicates the image of a corresponding function under the action of T _1 . We have 



max 



9 



< 



-i 



E\ and max 



h — Qn-l 



< 



1—1 1 



£2- From the d'Alembert formula by 



analogy with the standard proof of the stability of the Cauchy problem for the wave equation (see, 
e-g-, [SHI Sect. 4.3]) we obtain max^ t \ eA \u — un\ < 1 1 T 1 1| [s\ + e^b). Finally, (50) is obtained 

(u - 



from the following max^^ 



UN 



m ax( :Cit ) eA 



un)\ < ||T|| max (3 . ii)gi 



UN 



6 Approximation by the functions {(pk}kLo 

It follows from Proposition [28] that approximations of continuous functions by /- polynomials 



play a significant role in constructing approximate solutions of the Cauchy problem (43)-(44). 



The transmutation operator and the relation between the functions (p^ and the powers x made 



it possible to prove Proposition 27 showing that any continuous function may be approximated 
arbitrarily closely by finite linear combinations of the functions (fk- In this section we use the 
transmutation operators to extend some well-known results of approximation theory (see, e.g. jB], 
[7], [37]) onto approximations by the functions (ft and discuss different ways to construct such 
approximations for a given function. 

Denote by $ n , n = 0, 1, . . . the linear vector space spanned by the functions <pQ,...,ip n . It 



follows from Theorem 16 that the functions tpo,...,(p n are linearly independent, therefore the 
space <J> n is (n + l)-dimensional and the embedding <I> n C ^ n +i holds for any n. 
Define by 

Slid) = mm \\g - h n \\ 

the best approximation of a continuous function g by /-polynomials of degree n, i.e., by finite linear 



combinations Xlfc=o a k L Pk (Definition 10). Here || • || denotes the usual uniform norm on [—6, b]. Due 
to the embedding $ n C ^n+i the quantity £n{g) is monotone decreasing as n — > oo. Proposition 
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states that £n(g) — > 0, n — > oo for any function g G C[—b, b]. It is known in the approximation 
theory that for some functions g the convergence rate of £n(g) to zero may result to be arbitrarily 
slow. Nevertheless additional smoothness properties of the function g allow one to obtain more 
precise results on this convergence rate. 

Theorem 29 (Direct approximation theorem) Suppose the function g possesses on the seg- 
ment [—b, b] continuous f -derivatives of all orders up to the order k. Then for the best approximation 
by f -polynomials the following estimates hold for any n > k 

(— VllTHmaxlHT-illJTn 1 !!} 
gf( a ) < V 2 ^ " " U "'" "- y/ q l| 



and 



as n — > oo. 
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Proof. Consider the function g = T 1 g. As follows from Proposition 21 g G C fc [— 6, 6]. A variant 
of Jackson's theorem [H Chap. 4, Sec. 6] states that 

E ^ ± (n + l)n... 1 .(n-fc + 2) (y) fc|l ^ )|1 
where denotes the best approximation by algebraic polynomials of degree < n. Due to Propo- 



sition 
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we have ||fir^|| < max{||T -1 ||, 1 1 1 1 } • II^Csll- Now the first statement of the theorem 
follows from Theorem li~6l 

The second statement easily follows from another variant of Jackson's theorem [TJ VI. 2], |3T|, 
5.2.1]: if the function h € C k [— b, b], then for any n > k 

E n {h) < A(p n {x)) Lo(p n (x)), 



where p n (x) = V( fe sOfo+ jO _|_ u ^ := ^(/jWji) is the modulus of continuity of the derivative 
h^ k \ satisfying uj(t) — > 0, t — > 0, and the constant A does not depend on h and n. ■ 

Remark 30 A similar result holds under a weaker condition on the smoothness of the function g, 
namely, suppose that g possesses on the segment [—b, b] continuous f -derivatives of all orders up 
to the order k — 1 and the f -derivative of the order k — 1 is Lipschitz continuous on [—b,b], i.e., 
\d k _ l g(x) — d k _ l g{y)\ < M\x — y\ for some constant M and for every x,y G [—b,b]. Then there 
exists a constant C > such that 

C 

£l(g) < —u. for any n>k. 
n K 



The proof may be done similarly to the proof of Theorem 29 with the use of Jackson's theorem 
Chap. 4, Sec. 6] or \3% 5.2.4] an d the fact that if a function g is Lipschitz continuous, then the 
function g = T _1 g is Lipschitz continuous as well. 

The classical reasoning in the proof of an inverse theorem for the function g = T _1 g with the 
application of Markov's inequality and of an inequality for the derivative of the polynomial (see, 



e.g. |37| 4.8.7 and 6.2], [7, VII. 2]) allows us to prove a partial reverse statement of Theorem 29 
We show that the obtained convergence rate of the best approximations is close to optimal. 

Theorem 31 (Inverse approximation theorem) Suppose that the best approximations by f- 
polynomials of some function g satisfy for some integer number r and positive constants M and £ 
the inequality 

4(9) < Vn G N. (51) 

Then the function g possesses f -derivatives of order r in (—6, b) and f -derivatives of order at least 
[r/2] at the endpoints, where [■] denotes the integer part of a number. 



Proof. Consider the function g = T l g. As it follows from (51) and Theorem 16, there exists a 
sequence of polynomials P n such that 

M 

||?--Pn||<^ VneN, (52) 
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where M = M ||T 1 ||. Consider the series 

oo 

Pi(oo) + ^2( P 2 k +^ x )- p 2 k ( x ))- (53) 

k=0 

It is uniformly convergent due to the estimate 

M M 2M 

\P 2 k+i{x) -P 2k (x)\<\g- P 2 u\ + \P 2k+i -g\< -^y + 2(fc+1)(r+e) < ^iy, (54) 

and as it is easy to see, the sum of the series is equal to g. To finish the proof, we use two well-known 
inequalities for the derivative of the polynomial of order n defined on [—6, 6]. First, 

77 

\P' n {x)\ < j= = ||P n (x)[| 



and Markov's inequality 



yW — .)'- 



\P' n ix)\ < T \\P n (x)\\. 



From the first inequality and estimate (54) we obtain for any segment [— d, d] C (—6,6) 



, p (r) , v p (r), s| , 2MC d • 2( fc + 1 ) r 2 r + 1 MC d 

\p; k ' +1 (x) - p; k '{x)\ < 



2 fc+U-V ± 2 k \ > I — 2 fc ( r+£ ) 2 fc£ ' 

where the constant Cd depends only on r and the segment [— d, d]. The obtained estimate leads 



to the uniform convergence of the series of r-th derivatives of (53) and hence to the conclusion 
that g £ C r (—b,b). Similarly, the second inequality leads to the conclusion that g G C^ r ^ [— 6, 6]. 
Application of Proposition [21] finishes the proof. ■ 

Contrary to the L2" n orm, the problem of explicit finding of a polynomial of the best uniform 
approximation can be solved in some special cases only. But from a practical point of view the 
exact solution is not that necessary, it is enough to know a polynomial which is sufficiently close to 
the best one. Techniques such as least squares approximation or the Lagrange interpolation (with 
specially chosen nodes) work well though in general far from the best, see [33]. Below we briefly 
describe the iterative algorithm of E. Remez for constructing polynomials arbitrarily close to the 
best one. Even the zero step of the algorithm, the so-called Tchebyshev interpolation, usually gives 
better results then the Lagrange interpolation. For a detailed description of the algorithm with 
implementation details and all the required proofs we refer to [3T], [27], [3]. 

First we remind some definitions and statements related to Tchebyshev uniform approximations. 
See [B], [37] for details. 

A linear subspace V of C[—b, b] of (finite) dimension n + 1 is said to fulfill the Haar condition 
if it possesses the property that every function in V which is not identically zero vanishes at no 
more than n points of [—6, 6]. An equivalent condition is that the interpolation problem is uniquely 
solvable, i.e., for every set of n + 1 points x^ (k = 0, 1, . . . , n) in [—6, b] and every prescribed vector 
(z/0) 2/1) • • • j Dn) there exists a unique function h £ V such that 

h(x k ) = y k , k = Q,l,...,n. 



19 



If V is spanned by the functions ho,h\, . . . ,h n , another equivalent condition is that every determi- 
nant 



ho(xo) /ii(xo) ••• h n (x ) 
ho(x%) hi{xi) ... h n (xi) 



^0 



ho(x n ) hi(x n ) ... h n (x n ) 

for any distinct points xq, x%, . . . ,x n from [—b, b\. The Haar condition is necessary and sufficient 
for the unique solvability of the approximation problem. 

A system of linearly independent functions ho,h\, . . . ,h n is called a Tchebyshev system if 
the linear subspace spanned by these functions satisfies the Haar condition. 

Proposition 32 Let f be a real-valued non-vanishing continuous function on [—b,b]. Then the 



system of functions {(p^'kLo constructed by (14) is a Markov system, i.e., for any n £ No the first 
n+1 functions form a Tchebyshev system and the subspace <£ ra spanned by these functions satisfies 
the Haar condition. 

Proof. The proof by induction is straightforward using the fact that for the /-derivative the Rolle 
theorem holds. Also the result may be deduced from §3.11] if we observe that for the real- valued 
non-vanishing function / the system {ty?fc}]? =0 is a scaled Polya system. ■ 

Remark 33 As the following example shows, for f being a complex-valued function the Haar 
condition may fail for the subspaces <I> n . Consider f(x) = e lx . Then the first two functions ipt are 
cpo = e tx , <pi = since, and for large segments the function tp\ may have arbitrarily many zeroes. 

Assume that the function / and hence all functions (ph are real-valued (we briefly discuss the 
complex- valued case at the end of this section). 

The Remez algorithm is based on the Tchebyshev theorem with a generalization by de la Vallee 
Poussin which gives a characterization of the polynomial of the best approximation \37\ 2.7.3]. 

Theorem 34 (Tchebyshev's alternance theorem) If P n = Ylk=o c k^Pk is a polynomial with 
respect to some Tchebyshev system {(pk}2=o> 9 ^ a continuous function and Q is an arbitrary 
closed subset of the segment [a,b], then P n is the best approximation to g on Q if and only if the 
difference g{x) — P n (x) attains a maximum of its modulus on Q, with alternative signs, at least at 
n + 2 distinct points of the given set. 

Such set of n + 2 points is often called an alternant of the function g. An important consequence 
of this theorem is that for any continuous function g there exist exactly n + 2 distinct points 
£o, • • • , £n+i from [—6, b] such that the best approximation of g on the whole segment [—b, b] coincides 
with the best approximation on this so-called characteristic set of n + 2 points. The idea of the 
Remez algorithm is to construct iteratively subsets of [—b, b] each of them consisting of n + 2 points 
in such a way that on every step the value of the best approximation on the n + 2 points subset be 
increasing. 

In the case when the set Q consists of exactly n + 2 distinct points xo < X\ < . . . < x n+ \, the 
problem of determination of the best approximation polynomial P n of g on Q is exactly solvable 
and reduces to the solution of the system of n + 2 linear equations 

n 

Y J CkMx j ) + (-l) j E(g) = g(x j ), j = 0, 1, . . . ,n + 1 (55) 

k=0 
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for the coefficients Ck, k = 0, . . . , n and the value of the best approximation E(g) = £ n {9) on the 



set Q. The solution of the problem (55) for given points xq < x\ < . . . < x n+ \ and values of the 
function g(xj) in these points is also called Tchebyshev interpolation. Note that unlike the 
Lagrange interpolation, the resulted polynomial does not pass exactly through the given values of 
the function but the deviations of the polynomial from the given values are equal by absolute value 
at all points and differ only in sign. 

Let us describe the iterative algorithm of E. Remez. We are looking for an /-polynomial close 
to the one giving the best approximation £n(g) of a given function g by polynomials from *&n- 

We begin with a set Mq consisting of n + 2 distinct points —b < x^ < xf* < . . . < 



Corresponding to these points, using ( 55 ) we construct an /-polynomial of Tchebyshev interpolation 



go = Ylk=o c Wk ^ The function go(x) is the best approximation of g(x) on the set Mq. Denote 



the value E(g) obtained from (55) by Eo(g), and let Do := \\g — go\\- It follows from Theorem 34 
and from the observation that the best approximation on n + 2 points subset is not worse than the 
best approximation on the whole segment [—6, b] that 



\E {g)\ < £f(g) < D = \\g 



9o I 



Now either ||<?o — <?|| = Idols') I and we are done, or \\g — go 1 1 > 1-^0 (<?) I- The idea of E. Remez is 
to construct a new set M\ which again consists of n + 2 points, but for which the corresponding 
linear functional E\{g) has a larger magnitude than \Eo(g)\. 

There are two possibilities to define the set M\. The first is the so-called single exchange method. 
Exactly one of the points of Mo is replaced by a new point £ satisfying \g{Cj — <?o(£)l = \\9 ~ Soil- 
The point to be removed is chosen in such a way that the difference g — go alternates in sign at the 
points of the new sequence, it is not hard to derive an exact table of rules. Renumeration of the 
points according to their magnitudes produces the set M\. 

The second possibility is the general method of E. Remez. It involves simultaneous exchanges. 
The function ho '■= g — go possesses at least n + 1 zeroes zP\ k = 1,... ,n + l in the interval (—b,b) 



and 



~'k "* < z i+i < fc — 0, 1, ... ,n. 



Set = —b, £„°j^ 2 = b. Now in each interval = \z k , z^i] > k = 0, . . . , n + 1 we determine a 
ho{x^) > h (x) for all x G Jfc if sgnh (x^) = 1, 



point arjk such that 



and 

ho(x^) < h (x) for all x G J fc if sgnho(xy) = -1, 

that is, we are looking for a maximum if the difference between g and the previous approximation 
is positive, and for a minimum, if the difference is negative. Note that corresponding maxima and 
minima always exist. Here we assumed that Eo(g) ^ 0. If Eo{g) = the points x^) x are to be 
chosen as a sequence of points at which ho(x) has alternatively a maximum and a minimum. 

The iteration is repeated until the quantity — - — k - , characterizing the closeness of the 

Dk 

found /-polynomial to the best one, is not sufficiently small. 

Under the condition that in each of the sets M m+ i there is a point £ such that |/i m (£)l = ll^mll 
both the single and the general exchange algorithms converge to the best approximation. The 
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convergence speed is at least linear, i.e., there exists a constant q < 1 such that 



£l(g) - \E m+1 (g)\ < q(£f(g) - \E m (g)\) 

(see [27] for details) . Under some additional assumptions on the smoothness of the function g and 
the functions {(Pk}k=o an< ^ ^ ne number and type of extremal points of the difference h = g — g in 
[—6,6], where g is the /-polynomial of the best approximation, the convergence rate is quadratic 
[27\ Thm. 84]. I.e., for practical purposes only few iterations are required. 

As with any iterative algorithm, an important question is to choose properly a good initial set 
Mq. One of the possibilities is to consider the function g of the best least-square approximation 
to g with respect to the functions (po, . . . , <p n . It is known |27l p. 129] that if the difference g — g 
does not vanish identically on [—6, b] then it possesses at least n + 1 zeroes on [—6, 6], hence it has 
at least n + 2 alternating points of maxima and minima. The coordinates of these extremal points 
may be considered as the starting set Mq. 

Another possibility (see [27, 4.1 and 7.2]) is recommended if it is necessary to construct ap- 
proximations of several functions with respect to the same functions (po, . . . , <p n . We consider the 
problem of finding the best approximation ip of the function <p n +\ by the functions (po, . . . , ip n . The 
function s := <p n +i — <p is not identically zero and possesses exactly n + 2 extremal points. These 
extremal points form a good initial set for the Remez algorithm. In the case when the functions (p k 
coincide with the powers x k , the function s coincides (up to a constant factor) with the Tchebyshev 
polynomial T n+ i(x/b) and the extremal points are given by = —b cos k = 0, . . . , n + 1. 

It is worth mentioning that the approximation problem may be discretized and interpreted as a 
linear programming problem and solved by available software. We take a finite subset X C [—6, 6] 
consisting of points x\, . . . , xn, where N > n + 2. The condition 



max 



g{x) - ^2c k ip k (x) 



k=0 



E 



can be written as 



-E < g(xj) - ^2 cWkixj) < E, j = 1, . . . , N. 



k=0 



Our problem is to minimize the linear function l-E + 0-CQ + ... + 0-c n subject to 2N linear 
constraints 



3/1 



E + ^jC k ip{xj) > g 

k=0 
n 

E -^2c k <p(xj) > -g(xj 



k=0 



j = l,...,N 



1,...,N. 



The obtained problem can be solved by a variety of methods available for solving linear program- 
ming problems, see [32], [33] for details. 

At the end of this section we return to the case of the complex- valued function /. As was 



mentioned in Remark 33 the Haar condition may fail for the subspace <!>„. Even if the Haar 
condition holds, there is no immediate generalization of the Remez algorithm for the complex- 
valued case. The reason is that the Remez algorithm is based on the existence of a characteristic 
set of a function consisting of exactly n + 2 points. We remind that a subset X C [—6, b] is called 
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characteristic for the function g if the best approximation of g on the whole segment [—6, b] coincides 
with the best approximation on the subset X, but does not coincide on any proper subset of X. 
Contrary to the real-valued case in the complex-valued case a characteristic set may contain any 
number of points between n + 2 and 2n + 3, see [35]. There is no simple way to determine the 
number of characteristic points for a given function. What is more, the given function may have 
several characteristic sets containing different numbers of points. 

In the existing algorithms the discretized problem is considered and solved directly as a nonlinear 
optimization problem, e.g., a convex programming problem [I], [30], or the problem is transformed 
into a semi-infinite programming problem with the use of the fact that \h\ = max^ g [ 0j 27r) R- e ( e *^ ■ h) ■ 
The dual problem is considered and discretized for the second time with respect to the angle (j> 
and solved by the simplex method [JJ , [10] or by a Remez-like algorithm [12] , [13] , [36] , [9] . If the 
obtained approximation is not sufficiently close to the best one, the optimality criterium of the 
best approximation [33], |35j is reformulated as a system of nonlinear equations and the Newton 
iterations are used to improve the accuracy, see [10J , [36], [ID], [9] for details. 



7 Numerical examples 



In this section we present several numerical examples illustrating the application of the described 
results on generalized wave polynomials and approximation by functions {^fcjfcLo to numerical so- 



lution of the Cauchy problem ( 43 ) , ( 44 ) . On the first step the initial data g and h are approximated 
by /-polynomials and then the approximate solution of the problem (43), (44), the function un 



from Proposition 28 , is calculated on a mesh of points in the triangle from the figure in Section [5] 
and compared to a corresponding exact solution. All calculations were performed using Matlab in 
the machine precision. For the construction of the system of the functions ipk the following strategy 
was implemented using two Matlab routines from the Spline Toolbox: on each step the integrand 
is approximated by a spline using the command spapi and then it is integrated using f nint. This 
leads to a good accuracy, and the computation of the first 180-200 or even more functions ip^ 
proved to be a completely feasible task. In all the reported examples the number of subintervals in 
which the considered segment is divided when the integrand is approximated by a spline was 3000 
and the splines were of the forth order. In the presented numerical results we specify the parameter 
n which is the number of the calculated functions 



Example 35 Consider the Cauchy problem 

Du — c 2 u = 0, 



-b < x < b, t > 0, 



u(x, 0) = g(x) = cosh \J c 2 — Xfx, u t (x, 0) = h(x) 



The exact solution of this problem has the form 

1 



u(x, t) 



sin ct + cos Ait cosh \ 6 



\\x. 



c, Ai G C. 



(56) 
(57) 



The corresponding second-order ordinary differential equation (15), f" — c 2 f = admits a nonva- 
nishing solution f{x) = e cx , /(0) = 1. Based on this solution we construct n functions iff. defined 
by (14) and (ll)-(l3[! with xq = 0. The initial data for this example were chosen such that both g 



and h admit uniformly convergent on [—b,b] generalized Taylor series (see Subsection 3.2) whose 



23 



x 10 




1.2 

1 

0.8 
0.6 
0.4 
0.2 



x 10 



Figure 1: Graphs of \g — P 2 q\ (° n the left) and \h — Q\g\ (on the right) from Example 35 with the 
coefficients a& and fik obtained by (|48|) which in this case reduces to (58) and (59). 



expansion coefficients are known explicitly. Indeed, observe that g and h are solutions of the equa- 
tion v" — c 2 v = Xv with different values of the parameter X. In the case of g: A = —X\ and in 
the case of h: A = — c 2 . Since g(0) = 1 and g'(0) = due to Theorem^ the function g can be 
represented as follows 

g = Ul -cu 2 = ^ -pyr^* - (gfc + l)!^ 1 

where u\ and u 2 are defined by and the series are uniformly convergent on [—b, b] for any 
finite b. Thus, the coefficients olu from (45) have the form 

X 2k X 2k 
«2fc = (-l) fe ^jyy, «2fc+i = (~ 1 )* !+lc (2A; + 1)! ' fc = » 1 >— ( 58 ) 



Analogously we obtain 



~ (-c 2 ) k ^ (-c 2 ) h 

h= ul -c U2 =Y: w ^ - cx: ^2^^ 

k — k — 



and the coefficients (3^ from (45) have the form 

2k 2fc+l 

ft* = ( _1 )*(|j^|> /W - (-D^'^nji. *- 0,1.2.... (59) 

^4s an example let us take b = 2, c = 3 and Ai = 1. Consider the f -polynomials P n and Q n -i 
from Proposition 28 with n = 20 obtained by truncating the generalized Taylor series (45). Figure 
[7] depicts the distribution of the absolute error of such approximation of the functions g and h. One 
can observe that the apparently simplier function h = 1 is approximated much worse (10~ s against 
10~ 8 ) by the truncated generalized Taylor polynomial. 
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Figure 2: Graphs of \g — P\\\ (on the left) and \h — Qx$\ (on the right) from Example 35 with the 
coefficients otk and fik obtained by means of the Remez algorithm. 



The distribution of the absolute error of approximation of the solution of the Cauchy problem 
\5&) , (51) is presented on Figure^ Typically for an approximation based on a Taylor expansion 
(generalized or not) the absolute error increases with the distance from the center. 

Obviously, neither always the expansion coefficients of the generalized Taylor series of the initial 
data are available in a closed form nor always a continuous function is representable in the form of 
such a series. In Section^ several other possibilities for approximating functions by f -polynomials 
were discussed. In the present example alternatively to the generalized Taylor expansion we also 
apply the Remez algorithm (with the single exchange method). The developed computer program in 
Matlab establishes the corresponding value of n for approximating a function by an f -polynomial 
after which (for n + 1, n + 2, etc.) the approximation cannot be significantly improved limited 
by the machine precision. Thus, for the considered example the function g was approximated by 
Pn meanwhile h was approximated by Qi%. Figure^ depicts the distribution of the corresponding 
absolute error of approximation. The maximum value of the absolute error for g was of order 10 -9 
and for h - 10 -8 . 

The distribution of the absolute error of approximation of the solution of the Cauchy problem 
(56), (51) is presented on Figure^ The maximum absolute error of the approximate solution is 
of the order 10 -8 and to the difference from the solution computed previously with the use of the 
generalized Taylor coefficients here the distribution of the error over the domain is more uniform. 



Example 36 In this example we consider the same problem (56), (51), again with 6 = 2 and 
Ai = 1 but now with another value of c, c = 5i. Again we have the generalized Taylor coefficients 
in a closed form \5$) and (59) but application of the Remez algorithm encounters an obstacle. As 
the function f{x) = e cx is complex valued (see Remark 33) the same is true for the functions ipk 



meanwhile as was explained in Section [6] the Remez algorithm is directly applicable only to real 
valued functions. 

Here in order to find the coefficients of the f -polynomials from Proposition 28 by a distinct from 
the generalized Taylor formula method we solve the corresponding linear programming problem as 
explained at the end of Section^ 

How well the approximation based on the generalized Taylor formula works is shown in Figures 
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Figure 3: The distribution of the absolute error of the approximate solutions of the Cauchy problem 
from Example [35] computed according to Proposition 28 with the /-polynomials of order n = 20 



and generalized Taylor coefficients (on the left) and with the /-polynomials of order 11 (for g) and 
18 (for h) with the coefficients obtained by the Remez algorithm (on the right). 



[7] and [6] where the distribution of the absolute error of approximation is depicted for g, h and the 
solution u respectively in the case n = 50. 

The corresponding results obtained by solving the linear programming problem (again, due to 
the limitation of the machine precision we used n = 14 for the approximation of g and n = 25 for 
the approximation of h) are presented in Figures^ and^ 

Example 37 In our final example a variable coefficient equation is considered. Namely, we solve 
the following Cauchy problem 

Uu-x 2 u = Q, -b<x<b, t> 0, (60) 

u{x, 0) = g(x) = e x2 ' 2 (l + jT e" 5 '^ , (61) 

u t (x,0) = h(x) = xe x2/2 . (62) 

Its exact solution is given by the expression 

u(x,t) = g(x) cosh t + h(x) sinh y/3t 

The corresponding ordinary second-order equation has the form 

f" ~ x 2 f = 0. (63) 

The functions g and h are solutions of the equation 

v" — x 2 v = Xv 
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Figure 6: The distribution of the absolute error of the approximate solutions of the Cauchy problem 
from Example [36] computed according to Proposition 28 with the /-polynomials of order n = 50 



and generalized Taylor coefficients (on the left) and with the /-polynomials of order 14 (for g) and 
25 (for h) with the coefficients obtained by solving a linear programming problem (on the right). 



with A = 1 and A = 3, respectively and hence again we have in our disposal the possibility to write 



down the coefficients aj~ and (3k from (45) explicitly. We compute a particular solution f of (63) 



numerically using the SPPS method as described in (18], then compute the functions (fk and the 
corresponding coefficients ctk and (3k analogously to Example\35\ The result for n = 20 is presented 
by Figures^ and^ where the error of approximation of g, h and u is depicted. 

Application of the Remez algorithm delivers the following results. The functions g and h were 
approximated by f -polynomials of order 16 and 19 respectively and the distribution of the absolute 
error of the approximate solution of (60)— (62) is depicted on Figure [^j As can be observed with 
this relatively small number of functions ipk involved in the approximation a remarkable accuracy 
in the final solution is achieved of the order 10~ 14 . 
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